rm(list=ls())

pacman::p_load(pscl, tidyverse, haven)

# read in data
df <- read_csv('flagdata/flag_data.csv')

# model
summary(m1 <- zeroinfl(sale4x6 ~ pwht19zip_01 + 
                         pblk19zip_01 + pcollege19zip_01 + medinc19zip_01 + 
                         pprotectlf19zip_01 + totpop19zip_01 + popden19zip_01 + pct_r|
                         pprotectlf19zip_01 + totpop19zip_01, data=df))

# table
stargazer::stargazer(m1,
                     type='text',
                     covariate.labels = c(
                       'Pct White','Pct Black','Pct College','Median Income','Pct Protective Services',
                       'Total Population','Population Density', "Percept Republican Vote (16)"
                     ),dep.var.labels = c('Flag Sales'),
                     star.cutoffs = c(0.05,0.01,0.001), 
                     out =
                       'figure1B.html')

# hyp scenarios
xhyp <- seq(min(m1$model$pwht19zip_01),max(m1$model$pwht19zip_01),length.out=20)
yhatsave <- matrix(ncol = 20, nrow=1000)

# model matrix
newdat <- with(m1$model, 
                 data.frame(pwht19zip_01 = xhyp,
                            pblk19zip_01 = mean(pblk19zip_01),
                            pcollege19zip_01 = mean(pcollege19zip_01),
                            medinc19zip_01 = mean(medinc19zip_01),
                            pprotectlf19zip_01 = mean(pprotectlf19zip_01),
                            totpop19zip_01 = mean(totpop19zip_01),
                            popden19zip_01 = mean(popden19zip_01),
                            pct_r = mean(pct_r)))
coefs <- coef(m1)[1:9]
names(coefs) <- gsub('count_','',names(coefs))
vcovs <- vcov(m1)[1:9,1:9]
colnames(vcovs) <- gsub('count_','',colnames(vcovs))
rownames(vcovs) <- gsub('count_','',rownames(vcovs))

set.seed(12345)
simb <- MASS::mvrnorm(10000, coefs, vcovs)
mm <- cfMake(sale4x6 ~ pwht19zip_01 + 
         pblk19zip_01 + pcollege19zip_01 + medinc19zip_01 + 
         pprotectlf19zip_01 + totpop19zip_01 + popden19zip_01 + pct_r, data=df,nscen=20)
mm <- cfChange(mm, 'pwht19zip_01', x=xhyp, scen=1:20)
out <- loglinsimev(mm, simb)
out <- data.frame(x=xhyp, y=out$pe, l=out$lower, u=out$upper)
write_csv(out, 'pred_probs_flags_model.csv')

